function beta_spline=reg_demand_share_mean(user_input,w0_start_mat,set_constraint)

    %%%%%%%%%%%%%%%%%%%%%%%%%%
    % UNPACK USER
    %%%%%%%%%%%%%%%%%%%%%%%%%%
        
[~,~,~,~,~,n1,n2,k1,k2,~,~,p_grid_constraint,~]=user_unpack(user_input);

    %%%%%%%%%%%%%%%%%%%%%%%
    % INITIALIZE 
    %%%%%%%%%%%%%%%%%%%%%%%

sizew=(n1+k1)*(n2+k2);
grid_dim=size(p_grid_constraint,1);

    %%%%%%%%%%%%%%%%%%%%%%%
    % PREPARE NAG
    %%%%%%%%%%%%%%%%%%%%%%%

a=[];  

bl1w=-99*ones(sizew,1);
bu1w=99*ones(sizew,1);

if set_constraint==0,
    bl=bl1w;
    bu=bu1w;
end;

if set_constraint==1,
    bl_sl=-999*ones(grid_dim,1);
    bu_sl=0*ones(grid_dim,1);
    %
    bl=[bl1w; bl_sl];
    bu=[bu1w; bu_sl];    
end;

w0=w0_start_mat(:,1);
    
if set_constraint==0,
    ncnln=int64(0);
end;
if set_constraint==1,
    ncnln=int64(grid_dim);
end;

nclin=int64(0);
nn=int64(sizew);
cjac=zeros(max(1,ncnln),length(w0));
clamda=zeros(length(bl),1);
r=zeros(length(w0));
istate = zeros(size(bl,1), 1, 'int64');

[~,lwsav,iwsav,rwsav,ifail] = e04wb('e04uc');
[lwsav, iwsav, rwsav, ~] = e04ue('Print Level = 1', lwsav, iwsav,rwsav);
[lwsav, iwsav, rwsav, ~] = e04ue('nolist', lwsav, iwsav,rwsav);
[lwsav, iwsav, rwsav, ~] = e04ue('Major Print Level==0', lwsav, iwsav,rwsav);
[lwsav, iwsav, rwsav, ~] = e04ue('Col', lwsav, iwsav,rwsav);
[lwsav, iwsav, rwsav, ~] = e04ue('derivative level=3', lwsav, iwsav,rwsav);
[lwsav, iwsav, rwsav, ~] = e04ue('Hessian=yes', lwsav, iwsav,rwsav);
[lwsav, iwsav, rwsav, ~] = e04ue('Major Iteration Limit==6000', lwsav, iwsav,rwsav);
[lwsav, iwsav, rwsav, ~] = e04ue('Step Limit==0.1', lwsav, iwsav,rwsav);
[lwsav, iwsav, rwsav, ~] = e04ue('Minor Print Level==0', lwsav, iwsav,rwsav);
[lwsav, iwsav, rwsav, ~] = e04ue('Monitoring File==-1', lwsav, iwsav,rwsav);
[lwsav, iwsav, rwsav, ~] = e04ue('Verify Level == -1', lwsav, iwsav,rwsav);
[lwsav, iwsav, rwsav, ~] = e04ue('Stop Objective Check at Variable == 46', lwsav, iwsav,rwsav);
[lwsav, iwsav, rwsav, ~] = e04ue('Stop Constraint Check at Variable == 3', lwsav, iwsav,rwsav);
[lwsav, iwsav, rwsav, ~] = e04ue('Function Precision == 1e-12', lwsav, iwsav,rwsav);
  
    %%%%%%%%%%%%%%%%%%%%%%%
    % CALL NAG ROUTINE
    %%%%%%%%%%%%%%%%%%%%%%%

flag_continue=1;
ii=1;
ii_max=size(w0_start_mat,2);

while (ii<=ii_max) && (flag_continue==1),

        w0_start=w0_start_mat(:,ii);
        
        if set_constraint==0,
            [~, istate, ~, cjac, clamda, ~, ~, r, w, ~, lwsav, iwsav,rwsav, ~] = e04uc(a, bl, bu, '', 'objfun_mean', istate, cjac, clamda,r, w0_start, lwsav, iwsav, rwsav,'n', nn, 'nclin', nclin, 'ncnln',ncnln,'user',user_input);
        end;
        if set_constraint==1,
            [~, istate, ~, cjac, clamda, ~, ~, r, w, ~, lwsav, iwsav,rwsav, ~] = e04uc(a, bl, bu, 'confun_mean', 'objfun_mean', istate, cjac, clamda,r, w0_start, lwsav, iwsav, rwsav,'n', nn, 'nclin', nclin, 'ncnln',ncnln,'user',user_input);
        end;

        if ifail==0,
            flag_continue=0;
        end;
        ii=ii+1;
        
end;

    %%%%%%%%%%%%%%%%%%%%%%%
    % RESULTS
    %%%%%%%%%%%%%%%%%%%%%%%

beta_spline=NaN*w;

if ifail==0,
    beta_spline=w;
end;

    %
